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Gravitational waves at suitable frequencies can resonantly interact with a binary system, inducing 
changes to its orbit. A stochastic gravitational-wave background causes the orbital elements of the 
binary to execute a classic random walk, with the variance of orbital elements growing with time. 
The lack of such a random walk in binaries that have been monitored with high precision over long 
time-scales can thus be used to place an upper bound on the gravitational-wave background. Using 
periastron time data from the Hulse- Taylor binary pulsar spanning ~ 30 years, we obtain a bound of 
h c < 7.9 x 10~ 14 at ~ 10 -4 Hz, where h c is the strain amplitude per logarithmic frequency interval. 
Our constraint complements those from pulsar timing arrays, which probe much lower frequencies, 
and ground-based gravitational-wave observations, which probe much higher frequencies. Interesting 
sources in our frequency band, which overlaps the lower sensitive frequencies of proposed space- 
based observatories, include white-dwarf/supermassive black-hole binaries in the early/late stages 
of inspiral, and TeV scale preheating or phase transitions. The bound improves as (time span) -2 
and (sampling rate) _1//2 . The Hulse- Taylor constraint can be improved to ~ 3.8 x 10~ 15 with a 
suitable observational campaign over the next decade. Our approach can also be applied to other 
binaries, including (with suitable care) the Earth-Moon system, to obtain constraints at different 
frequencies. The observation of additional binary pulsars with the SKA could reach a sensitivity of 
fc c ~3x 10~ 17 . 

PACS numbers: 04.30.-w 04.30.Db, 95.30.Sf 



I. INTRODUCTION 

The first direct detection of gravitational waves (GWs) 
will be a landmark event. With the advent of the Ad- 
vanced Laser Interferometer Gravitational Wave Obser- 
vatory (LIGO) and Advanced Virgo network of detectors 
in ~ 2015, and the rapid progress of pulsar timing arrays 
(PTAs), it is likely that the next few years will see this 
breakthrough come to pass. 

In this paper, we will investigate an alternative ap- 
proach to GW detection, based on precision orbital mon- 
itoring of binary systems. The most promising binary 
systems are those with (at least) one pulsar member, 
although the way our method works and the frequen- 
cies that are probed differ significantly from the PTA 
method. PTAs probe gravitational waves that pass be- 
tween the pulsars and us, distorting the arrival times of 
what would otherwise be very regular pulses. We are in- 
terested instead in how background GWs interact with 
the orbital dynamics of a binary system. We emphasize 
that we are not so much interested in the emission of 
GWs by the binary - an important subject in its own 
right - as in the changes to its orbital parameters due to 
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scattering with some external GW background. 

The scattering is especially effective at GW frequen- 
cies that match a harmonic of the binary's orbital fre- 
quency, thereby inducing a resonance. For a circular 
orbit, resonance occurs at twice the orbital frequency. 
For eccentric orbits, resonance takes place to varying de- 
grees at all harmonics of the orbital frequency, starting 
from the fundamental. Note that this scattering process 
is quite special and not frequently discussed. A binary 
is inevitably losing energy over time by emitting grav- 
itational waves; the external GWs merely introduces a 
small modulation of this overall energy loss (for exter- 
nal GWs of a sufficiently small amplitude). Depending 
on the relative phase between the external gravitational 
wave and the binary, a constructive resonance slows down 
the energy loss, while a destructive resonance speeds it 
up. Immersed in a stochastic background of GWs, the 
binary orbital elements will thus execute a classic ran- 
dom walk on top of the secular decay due to its own 
GW emission. On average, the expected excursion of 
orbital elements due to GW scattering vanishes. How- 
ever, the variance of the excursion is non-zero and in fact 
grows over time, as in Brownian motion. Such orbital 
excursions can therefore be used to detect GWs, and the 
lack of such excursions can be used to place bounds. We 
again emphasize the fundamental difference between our 
approach and that of PTAs, since we are interested in 
measuring actual changes in the binary's orbital param- 
eters, rather than merely apparent changes due to the 
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presence of GWs along our line of sight. 

Versions of this idea have been discussed in pioneer- 
ing work by The subject lay dormant for many 
years, perhaps due in part to the focus shifting to PTAs 
in discussions of pulsars as a detection tool. We wish to 
revive the discussion by: (1) taking advantage of over 30 
years of precision monitoring of the famed Hulse- Taylor 
binary, recognizing that the rms orbital excursion due 
to GW scattering increases with time, (2) generalizing 
earlier work by computing the random walk of a binary 
orbit with arbitrary eccentricity, due to scattering by any 
stochastic GW background, and (3) finding the minimum 
variance estimator for the stochastic GW power spec- 
trum, given periastron time data; a main result is that 
the rms fluctuation in periastron time grows as (time 
span) 3 / 2 . The approach we propose can be thought of as 
an astronomical version of Weber's resonance bar Q. 

The only data we thoroughly analyze in this paper is 
from the Hulse-Taylor (HT) system PSR B1913+16. It 
provides a current constraint on the stochastic GW back- 
ground that is weaker than the one from the Dopplcr 
tracking of Cassini, and neither constraint is very restric- 
tive. However, future observations of similar binary sys- 
tems are expected to improve the constraint considerably. 
Our method can also be applied to other binary systems 
which have also been monitored with high precision over 
long time scales, such as the Earth-Moon and Earth-Sun 
systems. In Sec. V we discuss the prospects and chal- 
lenges of obtaining constraints from them. We will use 
the term "detector binary" as the generic descriptor for 
the systems of interest. 



II. FORMALISM 

We limit our focus to a GW background that is 
stochastic in nature, i.e. during the course of observation, 
the GW signal is not dominated by a single source with 
a definite phase, but rather arises from a multitude of 
sources, contributing to a signal that is statistically sta- 
tionary. For practical purposes, this means a Gaussian 
random field, although our calculation does not rely on 
Gaussianity. Because of our interest in PSR B1913+16, 
we are particularly interested in harmonics of its orbital 
frequency /ht = 3.6 x 1CP 5 Hz. 

This frequency corresponds to the Hubble scale at 
a temperature of about 100 GeV. Early universe pro- 
cesses, such as preheating after low scale inflation or 
bubble collisions at the electroweak phase transition, 
generate a stochastic GW background at these fre- 
quencies |^4llj. However, more promising sources of 
GWs at these frequencies may come from the later uni- 
verse, in particular from a large population of dou- 
ble white dwarf binaries in the early stages of inspi- 
ral, and from supermassive black-hole binaries in the 
late stages. These sources emit GWs at frequencies 
f r = 1.3 x 10~ 4 Hz [10 5 GMc- 2 /a} 3 / 2 [M Q /M], where f r 
is the frequency in the source's rest frame, a is the mean 



orbital separation, and M is the total mass of the binary. 
Estimates by jl2Hl6| suggest that the rms strain ampli- 
tude per logarithmic frequency from white dwarf binaries 
is roughly h c ~ 10" 20 - 10- 19 (/// HT ) _2/3 , where / is the 
orbital frequency of interest, and the number of sources 
within the frequency width of interest (see below) is suf- 
ficiently large to give a GW background that is Gaussian 
random to good approximation. 1 

Supermassive black-hole binaries constitute another 
promising source of GWs. However, the GWs generated 
by these systems at frequencies ~ 10~ 4 Hz would not be 
stochastic in character, since the number of black-hole 
binaries potentially resonating with HT is much smaller. 
The interaction of the detector binary with GWs from 
these sources is therefore better characterized by indi- 
vidual events. In the course of each event, the rela- 
tive phase between the source and the detector binary 
remains coherent; the detector orbital elements would 
therefore change in a secular rather than a stochastic 
fashion. Over time, the detector binary would encounter 
different, uncorrelated events, and thus there would still 
be a net random walk of sorts, if viewed over sufficiently 
long timescalcs. However, obtaining quantitative con- 
straints on GWs of such a character would require a dif- 
ferent calculation from the one presented here - a subject 
we hope to address in the future. In this paper, we fo- 
cus instead on the classic random walk effect, relevant 
for GWs from double white-dwarf binaries or the early 
universe. 

Let (X, Y, Z) define the frame of the detector binary, 
with the binary orbit lying in the X-Y plane. Let (x, y, z) 
define the frame of a particular gravitational wave train, 
with z being the incident direction. We can go from 
{X, Y, Z) to (x, y, z) by performing three consecutive 
Z — Y — Z Euler rotations. The last Euler rotation can 
be ignored since it is equivalent to rotating among the 
polarizations of the GWs, which we average over in any 
case. Thus, without loss of generality: 

x = (X cos 4> + Y sin cj>) cos 8 , (la) 
y = Y cos <f) — X sin <f) , (lb) 
z = (X cos <j) + Y sin </>) sin 9 , (lc) 

where <f> and 9 are two Euler angles. 

The induced relative acceleration A = A x x+A y y+A z z 
of a binary system due to a gravitational wave incident 



Gaussian randomness can be checked, for instance, by com- 
paring the connected 4th moment against the second moment 
squared. Demanding the former is small compared to the lat- 
ter is equivalent to requiring no. of sources X (A 2 ) 2 /(A 4 ) S> 1, 
where A is the amplitude of the gravitational wave from a given 
source. For our application, the number of sources ~ 10 while 
{A 2 ) 2 /{A i ) ~ 10" 5 . 
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from the z direction is given by 



1 /■;■ 



A x = -Rxoxox - RxOyoy = _- [h+x + h x y) , (2a) 



A y = -Ryoxox - Ryoyoy — 2 ( h x x - /i+y ) , (2b) 



A z = 0, 



(2c) 



where Rfi Va p is the Riemann tensor, and /i as the am- 
plitude of the gravitational wave strain, with + and x 
subscripts denoting the two polarization states. Eq. @ 
allows us to write down the time dependence of the en- 
ergy, 



dE 
~dt 



H (A x x + A y y + A z z) 



= g [ h +(xx - yy) + h x (xy + yx) ) , (3) 



where u = mi " 12 is the reduced mass of the binary. Note 
that this gives the energy change solely induced by the 
incoming gravitational wave, which is on top of its orig- 
inal energy loss due to emission. Analogous expressions 
can be written down for the angular and center-of-mass 
linear momenta (see Appendix). 

In Eq. (j3)), we can see that the quantities related to 
the binary motion are the quadrupole components. It is 
straightforward to obtain their Fourier expansions as 



X(t) 2 =X* + J2 Qxx(n) cos(27rn/t) , (4a) 

71=1 
OO 

Y(t) 2 =Y 2 + Y, Qyy{u) coB(27rn/t) , (4b) 

71=1 

OO 

X(t)Y(t) = X Y + Qxv(n) sm(2nnft) , (4c) 



where X a and Y a are constants, and / is the orbital fre- 
quency. The Q's are the quadrupole moments [l7j : 



Qxx(n) = — J n - 2 {ne) - 2eJ n ^ 1 (ne) 



+ 2eJ„+i(ne) - J rl+2 (ne) , 



4a 2 



(5a) 

Wv(n) = -Qxx(n) + J„(ne) , (5b) 

77/ 

Qxr(n) = - e 2 ^J„_ 2 ("-e) - 2J n (ne) + J n+2 (ne)^ 

(5c) 

where a is the semi-major axis, and e is the eccentricity 
(a = 1.95 x 10 6 km, e = 0.617 for PSR B1913+16). For 
our purpose, it is convenient to define the following 4 
quadrupole moments, which are more closely connected 



to the dynamics in the (x, y, z) frame: 

Qi(n) = (cos 2 6 cos 2 cp - sin 2 4>)Qxx(n) 

— (cos 2 <j) — cos 2 9 sin 2 4>)QYY(n) , (6a) 

Q 2 (n) = sin20(l + cos 2 8)Q XY (n) , (6b) 

Q3(n) = 2 cos cos 2^ Q_x-y(»i) , (6c) 



Q 4 (n) = cos#sin2(M <2xx(rc) - Qwin) j . (6d) 
For example, Eq. ([3]) can be rewritten as: 

f = ff>»/) (7) 

71=1 

■ Qi{n) h + sin(2irn ft) + Q2(n) h + cos(2irn ft) 



+ Qz{n) h x cos(27rn/i) + Qi{n) h x sin(27m/i) 



Note that ft, x and /i + have identical statistical proper- 
ties, and are uncorrelated. Strictly speaking, the orbital 
motion used on the right hand side of Eq. (J7J) to compute 
dE/dt should be the actual motion, accounting for both 
the orbital decay over time due to GW emission, and the 
orbital perturbation due to scattering with the external 
GWs. However, since both are very small effects, it is a 
very good approximation to use the unperturbed orbit. 

Depending on the phase of the incoming strain, dE/dt 
can take either sign. Averaging over an ensemble of 
stochastic GWs would yield a vanishing change in the 
orbital energy of a detector binary; to find an observable 
signature of GWs, we must therefore compute the energy 
variance. Let AE be the energy change over some period 
of time T. It can be shown that its variance takes the 
form: 



(™)2\ 



(8) 



i—l n—1 



where i labels the energy change associated with the 
quadrupole Qi. As an example, the 7 = 3 term is given 
by 



(AEi n)2 )= / dt I dt'Q 3 {nf(h x (i)h x {t')) 



T r T 



(9) 



o Jo 
cos(2irnft) cos(2Trnft') 

rj, 2 1 «6 h c (nf) 2 2 
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where we use: 



(Mt)MO) = <M*)MO) 
_i r df_ 2 , i2v f(t-t') 



(10) 



with h 2 , representing the (total) power spectrum per 
logarithmic frequency interval. We also assume T 3> 
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l/(27rn/), and use: 



dt e i27rft cos(27m/t) " « -$(/' - n/) . (11) 



A related delta function identity explains why different n 
modes do not mix in Eq. ((5]). The fact that different i's 
do not mix is partly due to the fact that the two different 
polarizations are uncorrelated, and partly due to the fact 
that the analog of Eq. (JTTJ for mixed sin and cos terms 
vanishes. Combining all terms, we have: 



(12) 



n=l 

Qi(n) 2 + Q 2 (n) 2 + Q 3 (n) 2 + Q 4 (n) 2 



Using the virial relation: E 



-[G 2 ( mi 



m 2 ) 2 / 2 7i" 2 /2] 1/3 M = ~2ir 2 f 2 a 2 n, and the period P 
1//, this can be rewritten as: 2 



(13) 



with 



, hc(nf) 2 
h c (2f) 2 



Qi{n) 2 + Q 2 {n) 2 + Q 3 {n) 2 + Q 4 {n) 2 (14) 



The Qi's depend on the incidence direction of the GWs, 
so we average over (6,<f>) to find the net effect. However, 
it is worth noting that, even without averaging, (AE 2 ) 
or (AP 2 ) vary by at most a factor of 2 across the sky. 

Eqs. (|12p and (|13[) make clear that only harmonics 
of the orbital frequency / contribute to the rms en- 
ergy/period change - the hallmark of a resonance effect. 
The singling out of these frequencies stems from delta 
functions like the one in Eq. (fTTj) . which has a width 
Af ~ 1/T, where T is the duration of integration. We 
are interested in T from weeks to years (3> the orbital 
period of 0.323 days for PSR B1913+16), corresponding 
to A/ - 1CT 9 - 1CT 6 Hz. GWs within this width of the 
harmonics would contribute to the random walk of the 
binary elements. 3 



2 Note that when the system is being perturbed by the external 
GWs, the virial relation does not strictly hold on an instan- 
taneous basis. However, when averaged over many orbits, we 
find the virial approximation of relating changes in energy to 
changes in period to be a very good one, with corrections sup- 
pressed by P/T; here the period is defined in an average sense 

p T -4- P ( T ~) 

i.e. 2-7T = J T uj(t)dt, as opposed to 2ir divided by the in- 

stantaneous angular velocity uj. 

3 The fact that the cumulative change in energy AE fluctuates, or 
random walks, can be understood as follows. As T varies, so does 
Af, which controls which sources of gravitational waves con- 
tribute to resonant scattering with the detector binary. Since the 
sources have uncorrelated phases, AE undergoes random kicks 
as the relevant source population varies. 



If the binary orbit is circular, only the n = 2 harmonic 
contributes, whereas for an eccentric orbit, all harmon- 
ics including n = 1 contribute in principle. In practice, 
the quadrupole moments Qi(n) decrease with n, and the 
expected rms strain h c drops with frequency, which coun- 
teracts the strong n 5 dependence in Eqs. (fT2|) and (fl"3|) . 
Assuming h c oc freq. -2 ^ 3 and the orbital parameters of 
PSR B1913+16, the dominant contributions come from 
n = {1, 4, 5, 3, 6, 7, 8, 2}, in order of importance, with the 
n = 1 mode contributing nearly as much as the other 
modes combined. Under the same assumptions, the di- 
mensionless amplitude A ~ 10. (Changing the spectrum 
to h c oc frcq.^ 1 would only change A to ~ 10.25.) 

Eq. $13} thus tells us that 



AR, 



10/i c (2/)i 



(15) 



which can be understood intuitively as follows. During 
each orbital period, the fractional change in period is 
roughly given by the strain h c . The cumulative rms 
change scales up by the square root of the number of 
periods y/T/P, as expected in a Brownian random walk. 
The extra factor of 10 depends on the details: the shape 
of the orbit and the spectrum of h c - we choose to nor- 
malize at 2/ (twice the orbital frequency) to facilitate 
comparison between different binaries, including circular 
ones. 



III. DATA ANALYSIS 

As a specific application of the generic method derived 
in the preceding section, we use PSR B1913+16 as the 
detector binary. The pulsar data carry a wealth of in- 
formation about the system. For simplicity, we focus on 
the periastron time, recognizing that stronger constraints 
could potentially be obtained by analyzing the full time- 
of-arrival data, which we leave for future work. The pe- 
riastron time data of PSR B1913+16 were published in 
[l8j . They consist of 27 periastron time measurements, 
spanning from 1974 to 2006, each of which was obtained 
from monitoring the system over approximately 2 weeks 
(and ~ 2 hours per day over those 2 weeks). Let us label 
these times Tj, with i ranging from 1 up to N (N = 27 
in our particular case). They can be modeled as follows: 



T l =f l + AT, 



(16) 



where f, is a smooth component, ATi is the excursion 
induced by GW scattering, and rii represents noise. The 
smooth component takes the following form: 



Ti — a + f3p{i) + jp(i)" 



(17) 



where p(i) tracks the number of periods, known to high 
accuracy; a is the zero-point (we choose p(l) — 0); /3 is 
essentially the period and would be exactly the period if 
there was not a small change over time; 7 quantifies the 



periastron shift due to the small change in the apparent 
period, induced by two smooth processes: (1) the famous 
decay of the orbit due to the emission of GWs, and (2) 
galactic acceleration of the system as a whole. Process 
(1) dominates over (2), though for our purpose there is 
no need to differentiate between them. Both are small 
compared to the zeroth-order effect i.e. yp(i)/j3 < 9 x 
10 -8 , and we will refer to j3 as the unperturbed period 
P; 7 can be thought of as PP/2 where P is the rate of 
change of the period due to (1) and (2). 4 

The fluctuations due to noise rii and due to GW scat- 
tering ATi arc uncorrclatcd, and their respective corre- 
lation matrices are 



(niUj) = CV 



(ATi AT j) = 



(18) 



where Cjj is the noise matrix, which we treat as diagonal 
using error bars from the data, and is defined as 



F 3 =^A 2 P 2 (min[p(t),p(j)}f 
(3max[p(i),p(j)) - mm[p(i) , p(j)} 



(19) 



To derive the expression for 
fact that AT, = 



(ATi AT j), we use the 
J^dtAP(t)/P, and compute 
(AP(t)AP(t')) using the same technique we used to com- 
pute (AP(t) 2 } in Eq. (13]). 5 We also take advantage 
of the useful fact that p(i) = J^T* dt/P(t), and the ex- 
pressions above can be derived by noting that AP, ATi 
and P are small quantities. Henceforth, to avoid clutter, 
we suppress the i,j indices and use bold-faced symbols 
to represent matrices or vectors, wherever no confusion 
would arise. 

It is worth pointing out that Eqs. (|T8|) and (fl~9| imply 
the rms fluctuation in periastron time, for i = j, is: 



AT rms ~6/ lc (2/ HT )P 



3/2 



(20) 



where we have abbreviated Ti as T, and P as P. This 
result can be roughly thought of as coming from scaling 
Eq. (|15[) up by T. In other words, the rms fractional 
change in period per period is roughly the strain h c ; after 
a number of periods given by T/P, the rms fractional 
change in period becomes ~ h c x y/T/P; since the peri- 
astron time is cumulatively dependent on t he per iod, the 
rms change in periastron time is ~ h c x ^JT/P x T. It 



In most of the paper, we simply use P to denote the unperturbed 
period, when there is no danger of confusion. 
We have implicitly assumed that the periastron time shift from 
external GWs is entirely due to their effect on the orbital period. 
In reality, there can be an (apparent) periastron time shift com- 
ing from center-of-mass linear momentum imparted by GWs, or 
from fluctuations in the orbital eccentricity. However, the former 
is suppressed by further powers of the orbital velocity, and the 
latter does not lead to a cumulative effect on the periastron time. 



is this rapid growth of AT rms with time span T that we 
exploit to obtain constraints on h c . 

The minimum variance estimator for T is [l9l ]: 



Test. =L(L T C- 1 L)- 1 L T C- 1 T 
where L is an TV x 3 matrix: 



1 p(l) p(l) 2 
1 p(2) p(2) 2 

1 p(N) p(N) 2 



(21) 



(22) 



The corresponding minimum variance estimator for the 
strain power spectrum h 2 at 2/ is 

h% est. = V (T T W T C- 1 FC- 1 WT - A) , (23) 

where the matrix W is an TV x TV matrix defined by 

W = 1 -L(L T C" 1 L)- 1 L T C^ 1 , (24) 

and rj and A are numbers defined by: 

77 = 1/Tr. [Cr 1 FCT 1 WF T W T ] , (25) 



A = Tr.[C- 1 FC- 1 WCW T ] . 



(26) 



In deriving the above expressions, we have assumed the 
variance of the estimators are dominated by the noise rii 
and not the signal ATj. It can further be shown that the 
variance of estimator h 2 cst is 

K cst .)-K est.) 2 = 2»? 



(27) 



Because of the subtraction of the noise power spec- 
trum (the A term), the estimated h 2 est can be neg- 
ative, though its ensemble average cannot. Note also 
the estimate hi est implicitly assumes the shape of the 
power spectrum, through the quantity A (see Eq. (fT4"l) ) 
in the definition of Fij. We will quote results assuming 
h c (f) cc /~ 2 / 3 . Assuming the power to go as / _1 would 
alter our results by a negligible amount compared to the 
uncertainties involved. 



IV. RESULTS 

From the periastron time data of PSR B1913+16 [lj], 
we find the power per logarithmic frequency interval at 
2/ht = 7.2 x 1(T 5 Hz to be h 2 c cst = -6 + 7.4 x 1CT 27 . 
This is derived from fluctuations in the data after fitting 
and removing a smooth quadratic (see Eq. (|17jl ). From 
this, we derive a 95% upper limit of h c < 7.9 x 10~ 14 at 
7.2 x 1(T 5 Hz. 6 



6 This upper limit corresponds to the value of h c such that the 
probability of observing h^. cst at —6 X or less is 5%. 
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Pulsars are known to have glitches, which are recogniz- 
able by abrupt changes in the spin period. The periastron 
timing data exhibit a large excursion around the glitch of 
May 2003, which we have removed from the above anal- 
ysis. Had we included that data point, the results would 
not have changed significantly, as we would instead find 
K est . = -3 + 7.4xl0- 27 . 



V. DISCUSSION 

Our GW constraint from the random walk of binary or- 
bital elements raises several interesting issues. First, the 
constraint is a conservative one; namely, it is based upon 
a bound on fluctuations of the periastron data around a 
smooth curve. If there are additional sources of fluctu- 
ations other than scattering from the GW background, 
accounting for them would only strengthen our bound. 
However, a more thorough understanding of the possi- 
ble sources of fluctuations would be necessary if one were 
to claim a detection of the GW background. Possibili- 
ties include glitches and tidal effects. Glitches are dis- 
tinguished by accompanying fluctuations in the pulsar 
spin period, which are unlikely to have been caused by 
stochastic GWs. Tides cause secular changes of the bi- 
nary dynamics rather than stochastic changes. To isolate 
the GW signal from other random or near-random pro- 
cesses, one can also take advantage of the well defined 
shape predicted for the two-point function of the perias- 
tron time fluctuations (see Eq. (|18|l). 

The only other direct bound on GWs at a frequency 
~ 10~ 4 Hz comes from Doppler tracking of the Cassini 
spacecraft [20"j , which is roughly an order of magnitude 
more stringent than our constraint. 7 Our bound is rather 
weak, especially when compared to the expected GW 
background from white-dwarf binaries, as shown in Fig . 
1. The expected GW background is taken from [l4| . 
which is consistent with estimates such as [l5l . [l6l ]. though 
with a large uncertainty. Also shown is the expected sen- 
sitivity for the proposed eLISA/NGO/SGO detector dJ. 

Thus an important question is: how much do we expect 
the GW bound to improve from future observations of 
PSR B1913+16 and other binaries? To guide our think- 



In our method, most of the constraining power of the data 
comes from the data points that are furthest apart i.e. T ~ 32 
years, which corresponds to a fairly narrow frequency window of 
A/ ~ 10~ 9 Hz. If one were to extrapolate this amplitude of he 
(7.9 X 10 — 14 ) to a smooth spectrum over a broad bandwidth, one 
would obtain Qgw > 1 which we know is ruled out by cosmolog- 
ical observations already. From this point of view, our bound is 
certainly weak. But it should be kept in mind that our bound on 
h c applies strictly within a narrow frequency window, for which 
there is no useful cosmological bound. Note also that our bound 
is weak compared to the bound from big bang nucleosynthesis 
(BBN). However, once again, the BBN bound assumes a broad 
(scale invariant) spectrum of GWs, and it bounds GWs in the 
early universe and not from late time astrophysical sources. 



ing, we observe that the sensitivity of our method to the 
GW background scales as 

h c ~ 5 (^) N-W (^f-y 3 ^ (28) 

fST\ 1/2 {T tot .y 2 

where ST is the accuracy of each periastron time mea- 
surement, N is the number of such data points, P is the 
orbital period, T tot . is the total time span, and np is the 
number of periods between consecutive periastron time 
measurements, so that rip 1 is the sampling rate. 

The minimalist approach would be to simply lengthen 
T to t.- Assuming the same rate of sampling as before (one 
periastron time data point per year), out to the year 
2022, would push the sensitivity on h c to 8.9 x 10~ 15 , 
comparable to the Cassini bound. This assumes ST ~ 
3 x 10 -8 day, which is about the level of accuracy towards 
the later years of the periastron data we analyzed [l8| . 




frequency [Hz] 



FIG. 1: Current constraint and future sensitivity on h c , the 
square root of the strain power spectrum per logarithmic fre- 
quency. The uppermost line with the (red) square marker 
shows the current constraint from analyzing ~ 30 years of pe- 
riastron time data from PSR 1913+16. The middle line with 
the (green) circle shows the expected sensitivity from the same 
system, if it continues to be observed ~ 2 hours per day for 
~ 2 weeks each year until the year 2022. The bottom (blue) 
line shows the sensitivity from a hypothetical observational 
campaign employing the SKA to observe 100 binary systems 
(see text for details). In each case, the solid line going through 
the dot represents a (freq.) _2//3 spectrum which is assumed in 
deriving the bound. For comparison, we also show the ex- 
pected strain sensitivity for eLISA/NGO/SGO [1| (dashed 
line), and the expected signal strength from white dwarf bi- 
naries (dotted line) [13 ] . We caution that the white dwarf 
binary background has large uncertainties. 

A more ambitious approach would be to increase the 
sampling as well. Recall that the data we analyzed came 
from ~ 2 weeks of observations per year during which 
the pulsar was observed for only ~ 2 hours per day. This 
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is a fairly sparse sampling. What if we increase the sam- 
pling to one periastron time data point every 2 weeks 
(from 2012 to 2022)? The sensitivity on h c then becomes 
3.8 x 10~ 15 . Increasing the sampling is not as effective 
as increasing the total time span, but there is still some 
useful gain. 

To go beyond this, let us consider the possibilities of- 
fered by the Square Kilometer Array (SKA). The SKA 
is expected to find hundreds of binaries with a pulsar 
member [22|]. Let us assume 100 binaries, with an or- 
bital period of around 0.1 day. 8 The SKA also has a 
higher sensitivity than existing instruments. How much 
this translates into an improvement in the pulsar timing 
residual depends on how important the pulse jitter is, but 
an order of magnitude or so improvement is conceivable 
[23|. 9 Let us assume ST ~ 10 -9 day. (As an exam- 
ple, the periastron time data had improved in accuracy 
by almost two orders of magnitude from 1974 to 2006.) 
For the sampling, let us use N = 100 data points in the 
time span of 15 years. The projected sensitivity on h c 
becomes ~3x 10~ 17 . This is also shown in Fig. 1. 

There are a few open questions that remain to be ex- 
plored. One is whether an even stronger GW bound 
can be obtained by analyzing the time-of-arrival data di- 
rectly, as opposed to the periastron time data. There is 
a wealth of information in the time-of-arrival data, only 
a small fraction of which is captured by the periastron 
time. The question is whether, as far as the impact of ex- 
ternal GWs on the test binary is concerned, most of the 
information is already contained in the periastron time 
data. For instance, we have not used any information 
about changes to the orbital eccentricity due to scatter- 
ing with the external GWs, which can be deduced from 
the changes in energy and angular momentum (see Ap- 
pendix). How much can our constraint improve if we use 
such information as well? Another interesting question 
is what current level of constraint we can obtain from 
other binary systems. Two in particular come to mind: 
the double pulsars PSR J0737-3039A/B, and the Earth- 
Moon system. A rough estimate for the double pulsar 
system, discovered in 2003, can be obtained by using 
ST ~ 2 x 10~ 7 day, P ~ 0.1 day, T tot . ~ 9 years, N ~ 30 
p4| . giving a current sensitivity to h c of ~ 3x 10~ 13 . This 
bound will improve more rapidly than the bound from 
the HT binary, since it has been observed for a shorter 
time. The Earth-Moon system is sensitive to the GW 



We envision binaries spanning a range of periods. Each thus 
probe the GW background at a different frequency. Assuming 
a spectrum for the background, we can bound a single number, 
i.e. the amplitude of the spectrum, using all the binaries. We 
assume the external GWs at the 100 binaries can be treated as 
uncorrclatcd. If there is some overlap in frequencies within the 
relevant resonant widths, there is the interesting possibility of 
cross-correlating excursions between binaries, which we leave for 
future work. 

In cases where pulse jitter is important, improvement can only 
be achieved by longer integration. 



background at a very different frequency: ~ 8x 10~ 7 Hz. 
Laser-ranging to the moon can measure the Earth-Moon 
distance down to ~ 15 mm, which is a fractional accuracy 
of about 4 x 10 -11 [25j . The main hurdle to obtaining 
accurate constraints on GWs is the need to model many 
geophysical effects of both the Earth and the Moon. A 
conservative bound can be obtained as long as one does 
not over-fit the lunar-ranging data. 

Let us close by noting that our calculation applies 
strictly to a stochastic background. The case of super- 
massive black hole mergers needs to be separately con- 
sidered, since their relative scarcity implies a GW signal 
more in the form of individual events, each of which is 
coherent. We hope to explore this in a future paper. 
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Appendix A: Angular momentum 

For completeness, we provide equations that describe 
the change of angular momentum due to scattering with 
the external GWs. 

^tt = M (~ za y) = \l l (h+yz - h x xzj , 

^ = n (za x ) = ^fj, (h+xz + h x yzj , 
dL z 

—jj- = [i (xa y - ya x ) 

= ^i(-2h + xy + h x (x 2 -y 2 )'} , (Al) 
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Applying the same procedure as we did for the energy, The variance in the magnitude of L is given by 
the variance in the change in angular momentum is: 

4 oo 

((AL X ) 2 ) = {(ALy) 2 ) = yT/V5> 3 M»/) 2 

n— 1 

(A2a) 

Q 5 (n) 2 + Q 6 (n) 2 + Q 7 (n) 2 + Q 8 (n) 2 ) , ^ ^ _ ^ 

((AL Z ) 2 ) = ^7r 4 r/V E ^ 3/i c(«/) 2 (A2b) = |£r 2 ((i*AL a + LyALy + L Z AL Z ) 2 ) 



2 

n=l 

2 , /-> /„\2 , /„\2 , / \2 



Qi(n) 2 + Q 2 (n)^ + Q 3 (n)' + Q*(n) 
where Q$,Q§,Q7, Qs are defined by 



^-^((AL^ + L^AL,,) 2 ) 
+ L 2 <(AL S ) 2 )) 

= sin 2 9((AL X ) 2 ) + cos 2 6>((AL Z ) 2 > . (A4) 



Qs(^i) = sin 20 ^Qxx(n) cos + Qyy(n) sin 0j , 

Qe(n) — sin 20 sin 26* QxY( n ) , (A3a) 

Q 7 (n) = 2 sin 6 cos 2^ QxyO) , (A3b) 

Q 8 (n) = sin sin 20 [ Q X x (n) - Qyy (n) j . ( A3c) 
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